/* -*- mode: c++ coding: utf-8 -*-

  This file is part of the Life library

  Author(s): Christophe Prud'homme <christophe.prudhomme@ujf-grenoble.fr>
       Date: 2008-02-07

  Copyright (C) 2008-2009 Universite Joseph Fourier (Grenoble I)

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
/**
   \file code1d.cpp
   \author Christophe Prud'homme <christophe.prudhomme@ujf-grenoble.fr>
   \date 2008-02-07
 */
/** include predefined life command line options */
#include <life/options.hpp>

/** include linear algebra backend */
#include <life/lifealg/backend.hpp>

/** include function space class */
#include <life/lifediscr/functionspace.hpp>

/** include helper function to define \f$P_0\f$ functions associated with regions  */
#include <life/lifediscr/region.hpp>

/** include integration methods */
#include <life/lifepoly/im.hpp>
/** include gmsh mesh importer */
#include <life/lifefilters/gmsh.hpp>

/** include exporter factory class */
#include <life/lifefilters/exporter.hpp>

/** include gmsh generator for tensorized domains */
#include <life/lifefilters/gmshtensorizeddomain.hpp>

/** include  polynomialset header */
#include <life/lifepoly/polynomialset.hpp>

/** include  the header for the variational formulation language (vf) aka FEEL++ */
#include <life/lifevf/vf.hpp>

#include <fstream> 

/** use Life namespace */
using namespace Life;
using namespace Life::vf;

/**
 * This routine returns the list of options using the
 * boost::program_options library. The data returned is typically used
 * as an argument of a Life::Application subclass.
 *
 * \return the list of options
 */
inline
po::options_description
makeOptions()
{
    po::options_description code1doptions("Code1D options");
    code1doptions.add_options()
        ("hsize", po::value<double>()->default_value( 0.5 ), "mesh size")
        ;
    return code1doptions.add( Life::life_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Life::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
    AboutData about( "code1d" ,
                     "code1d" ,
                     "0.2",
                     "Code 1D for -Delta u = f",
                     Life::AboutData::License_GPL,
                     "Copyright (c) 2010 Universite Joseph Fourier");

    about.addAuthor("Christophe Prud'homme", "developer", "christophe.prudhomme@ujf-grenoble.fr", "");
    return about;

}


/**
 * \class Code1D
 *
 * Code1D Solver using continuous approximation spaces
 * solve \f$ -\Delta u = f\f$ on \f$\Omega\f$ and \f$u= g\f$ on \f$\Gamma\f$
 *
 * \tparam Dim the geometric dimension of the problem (e.g. Dim=1, 2 or 3)
 */
template<uint N>
class Code1D
    :
    public Application
{
    typedef Application super;
public:

    //! Polynomial order \f$P_2\f$
//    static const uint16_type Order = 5;
    ///static const uint Order = N;

    //! numerical type is double
    typedef double value_type;

    //! linear algebra backend factory
    typedef Backend<value_type> backend_type;
    //! linear algebra backend factory shared_ptr<> type
    typedef boost::shared_ptr<backend_type> backend_ptrtype;


    //! sparse matrix type associated with backend
    typedef backend_type::sparse_matrix_type sparse_matrix_type;
    //! sparse matrix type associated with backend (shared_ptr<> type)
    typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    //! vector type associated with backend
    typedef backend_type::vector_type vector_type;
    //! vector type associated with backend (shared_ptr<> type)
    typedef backend_type::vector_ptrtype vector_ptrtype;

    //! geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1
    typedef Simplex<1> convex_type;
    //! mesh type
    typedef Mesh<convex_type> mesh_type;
    //! mesh shared_ptr<> type
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    //! the basis type of our approximation space
    //typedef bases<OrthonormalPolynomialSet<Order,Scalar> > basis_type;
    typedef bases<Lagrange</*Order*/N,Scalar,PointSetFekete> > basis_type;

    //! the approximation function space type
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    //! the approximation function space type (shared_ptr<> type)
    typedef boost::shared_ptr<space_type> space_ptrtype;
    //! an element type of the approximation function space
    typedef typename space_type::element_type element_type;

    //! the exporter factory type
    typedef Exporter<mesh_type> export_type;
    //! the exporter factory (shared_ptr<> type)
    typedef boost::shared_ptr<export_type> export_ptrtype;

    /**
     * Constructor
     */
    Code1D( int argc, char** argv, AboutData const& ad,
               po::options_description const& od )
        :
        super( argc, argv, ad, od ),
        M_backend( backend_type::build( this->vm() ) ),
        //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
        meshSize( this->vm()["hsize"].template as<double>() ),
        exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
    {
    }

    /**
     * create the mesh using mesh size \c meshSize
     *
     * \param meshSize the mesh characteristic size
     * \return a mesh of an hyper cube of dimension Dim
     */
    mesh_ptrtype createMesh( double meshSize );

    /**
     * run the convergence test
     */
    void run();

private:

    /**
     * solve the system \f$D u = F\f$
     *
     * \param in D sparse matrix
     * \param inout u solution of the system
     * \param in F vector representing the right hand side of the system
     */
    void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );


    /**
     * export results to ensight format (enabled by  --export cmd line options)
     *
     * \param u the function to save in ensight format
     */
    void exportResults( element_type& u );

private:

    //! linear algebra backend
    backend_ptrtype M_backend;

    //! mesh characteristic size
    double meshSize;

    //! exporter factory
    export_ptrtype exporter;
}; // Code1D

template<uint N>
typename Code1D<N>::mesh_ptrtype
Code1D<N>::createMesh( double meshSize )
{
    /** instantiate a new mesh */
    /** \code */
    mesh_ptrtype mesh( new mesh_type );
    /** \endcode */

    //! generate a tensorized domain (hyper cube of dimension Dim)
    /** \code */
    GmshTensorizedDomain<convex_type::nDim,
                         convex_type::nOrder,
                         convex_type::nDim,
                         Simplex> td;
    td.setCharacteristicLength( meshSize );
    td.setX( std::make_pair( 0, 2 ) );
    std::string fname = td.generate( convex_type::name().c_str() );
    /** \endcode */

    //! importer the mesh generated by td
    /** \code */
    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );
    /** \endcode */

    return mesh;
} // Code1D::createMesh


template<uint N>
void
Code1D<N>::run()
{

    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }


    int MyOrder=N;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % MyOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    mesh_ptrtype mesh = createMesh( meshSize );


    //The function space and some associated elements(functions) are then defined
    space_ptrtype Xh = space_type::New( mesh );
    element_type u( Xh, "u" );
    element_type v( Xh, "v" );
    
    value_type pi = M_PI;
    AUTO( g, exp(-Px()) );
    AUTO( f, -exp(-Px()) );
    
    vector_ptrtype F( M_backend->newVector( Xh ) );

    form1( _test=Xh, _vector=F, _init=true ) =
        integrate( elements(mesh), _Q</*Order*/N+5>(),
                   f*id(v) );
    
    form1( _test=Xh, _vector=F ) +=
      integrate( markedfaces(mesh,1), _Q<2*N /*Order*/>(),id(v) );

    F->close();



    sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );

    form2( Xh, Xh, D, _init=true ) =
        integrate( elements(mesh), _Q<2*N /*Order*/>(),
                   gradt(u)*trans(grad(v)) );

    // weak dirichlet conditions treatment for the boundaries marked 1 and 3
    // -# assemble \f$\int_{\partial \Omega} -\nabla u \cdot \mathbf{n} v\f$
    // -# assemble \f$\int_{\partial \Omega} -\nabla v \cdot \mathbf{n} u\f$
    // -# assemble \f$\int_{\partial \Omega} \frac{\gamma}{h} u v\f$

    form2( Xh, Xh, D ) +=
      integrate( markedfaces(mesh,2), _Q<2*N /*Order*/>(),
                 idt(u)*id(v) );
    D->close();

    this->solve( D, u, F );

//    std::cout<<"\nnDof="<<u.nDof()<<"\n";
//    std::cout<<"\nnElem="<<u.functionSpace()->mesh()->numElements()<<"\n";
      
//    for (uint i=0;i<u.nDof();++i) {
//      std::cout<<" "<<u(i); }

/*
    std::cout<<"\nles dofs=\n";
    for (uint e=0;e<u.functionSpace()->mesh()->numElements();++e)
      for (uint i=0;i<u.nDof();++i) {
//        std::cout<<" "<<u.functionSpace()->dof()->localToGlobal( e, i );
        u.functionSpace()->mesh()->point(i);//.point(0).showMe();
    }
*/

    //-----------------------------------------------------------------------------------//

    std::ostringstream Nstring;Nstring << N;
  
    std::string name("sol1d_N"+Nstring.str()+".txt" );

    std::ofstream fichier;
    fichier.open(name.c_str(), std::ios::out);

    typename space_type::dof_type::dof_points_const_iterator it_dofpt = u.functionSpace()->dof()->dofPointBegin();
    typename space_type::dof_type::dof_points_const_iterator en_dofpt = u.functionSpace()->dof()->dofPointEnd();

    for( ; it_dofpt != en_dofpt; ++it_dofpt  ) 
      fichier << boost::get<0>( *it_dofpt )(0) << " "<< u(boost::get<1>( *it_dofpt ))<<"\n";

    fichier.close();

    //-----------------------------------------------------------------------------------//

    //! compute the \f$L_2$ norm of the error
    double L2error2 =integrate(elements(mesh), _Q<2*N /*Order*/>(),
                               (idv(u)-g)*(idv(u)-g) ).evaluate()(0,0);
    double L2error = math::sqrt( L2error2 );

    this->changeRepository( boost::format( "%1%/%2%/Erreur/h_%3%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % this->vm()["hsize"].template as<double>()
                            );

    std::ofstream fichierErr;
    fichierErr.open("Err.txt", std::ios::app);

    fichierErr<< N << " " << L2error << "\n" ;

    fichierErr.close();

//    std::cout<<"\n\n||error||_L2=" << L2error << "\n";
//    Log() << "||error||_L2=" << L2error << "\n";

    //-----------------------------------------------------------------------------------//

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % MyOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    this->exportResults( u );

} // Code1D::run

template<uint N>
void
Code1D<N>::solve( sparse_matrix_ptrtype& D,
                  element_type& u,
                  vector_ptrtype& F )
{
    //! solve the system, first create a vector U of the same size as
    //! u, then call solve,
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    //! call solve, the second D is the matrix which will be used to
    //! create the preconditionner
    M_backend->solve( D, D, U, F );
    //! copy U in u
    u = *U;
} // Code1D::solve


template<uint N>
void
Code1D<N>::exportResults( element_type& U )
{
    Log() << "exportResults starts\n";

    exporter->step(0)->setMesh( U.functionSpace()->mesh() );

    exporter->step(0)->add( "u", U );

    exporter->save();
} // Code1D::export



/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
    Code1D<1> code1dN1( argc, argv, makeAbout(), makeOptions() );
    Code1D<2> code1dN2( argc, argv, makeAbout(), makeOptions() );
/*    Code1D<3> code1dN3( argc, argv, makeAbout(), makeOptions() );
    Code1D<4> code1dN4( argc, argv, makeAbout(), makeOptions() );
    Code1D<5> code1dN5( argc, argv, makeAbout(), makeOptions() );
    Code1D<6> code1dN6( argc, argv, makeAbout(), makeOptions() );
    Code1D<7> code1dN7( argc, argv, makeAbout(), makeOptions() );
    Code1D<8> code1dN8( argc, argv, makeAbout(), makeOptions() );
    Code1D<9> code1dN9( argc, argv, makeAbout(), makeOptions() );
    Code1D<10> code1dN10( argc, argv, makeAbout(), makeOptions() );  
*/
    code1dN1.run();
    code1dN2.run();/*
    code1dN3.run();
    code1dN4.run();
    code1dN5.run();
    code1dN6.run();
    code1dN7.run();
    code1dN8.run();
    code1dN9.run();
    code1dN10.run();
    */
}
